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GENERALIZED DUFORT-FRANKEL SPECTRAL METHODS 


By 

Liviu Lustman* 

ABSTRACT 

This report presents an explicit time-advancing scheme for the 
spectral solution of parabolic equations. Several two-dimensional 
examples are considered, including convection-diffusion and nonlinear 
problems, under various boundary conditions. Numerical evidence demon- 
strates the efficiency and accuracy of the spectral approach. 

This work was performed in collaboration with David Gottlieb, Tel Aviv 
University/ Institute for Computer Applications in Science and Engineering. 

GENERALIZED DUFORT-FRANKEL SCHEMES 
Introduction 

Gottlieb and Gustafsson (ref. 1) presented an explicit, unconditionally 
stable formula for time-advancing parabolic equations, while preserving 
high accuracy in space. The present work treats, as a special case, 
the Tschebyscheff collocation spectral method as the spatial derivative 
approximation. Large At may be employed, much above the limit OCN"**) 
needed for straightforward explicit methods which use Tq, Ti . . . Tj^. 

Thus the computation becomes efficient, while spectral accuracy is still 
maintained. Let us review the generalized Dufort-Frankel method for the 
simplest diffusion equation: 


The original Dufort-Frankel scheme is 
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Old Dominion University, Norfolk, VA 23508 



n^l 
ir 


2Lt 


s a- 


a 


( 2 ) 


n-1 


u. 


n 


♦ u 


n 


ti 


o n 
- 2u. 


J. . 




(AxV 


, n+1 

'“i 


n-1 
♦ u. 

3 




As usual, u? “ u(jAx, nAt) . The first term on the right-hand side is the 
simplest approximation to u^(jAx, nAt). It may be replaced by other 
approximations, e.g. higher order finite difference, finite elenent, or 
spectral. The second term on the right-hand side is also modified by a 
multiplier y > 0« Large enough values of y ensure unconditional stability 
(as proved for various cases in ref. 1). 

Finally, the method we shall employ is 
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Here u = u(x., nAt), where x.( 0 < j < N) is the Tschebyscheff 
J J J ~ ~ T, 

collocation mode; Ax may be taken to be 1 - cos — (minimal distance 

• ^ 

between nodes) or, for better accuracy: 


Ax = Ax. * x. 
J 3 
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Obviously (u^^)j 
explicitly for u' 


is computed spectrally, and equation (3) is solved 
n+1 

j 


It is important to mention that equation (3) is consistent, within 
O(At^), with a wave equation: 
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where the speed c is given by 



However, equations (1) and (5) share the same steady states as t -*• <»; 
also, if At = O(Ax^), equation (3) becomes consistent with equation (1) 
within O(At'). 
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The generalization to two dimensions: 
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It is also possible to replace AXj and by their miniaial values. 

Note the simplicity in the treatment of the mixed term u 

xy 

It is convenient to define ••fluxes” for equation (5) since these 
are the quantities actually treated by the numerical scheme: 

3 ^ 3 

t ax 3y ^ 


( 8 ) 


( 9 ) 


with fluxes f, g: 


f * j 

g = <^3Uy ♦ 7 ° 2 \ 

Note that there is nonuniqueness in the definition of f and g: 
for example, one might take 

t . o,u^ 

g « 02U^ + 03 Uy 

and still end up with equation (7). 

The solution is implemented as follows: 

Algorithm 

Cl) The first two levels u9j^, ujj^ are stored (u° is, of course, 
the initial data, but u^ must be obtained by some different method). 
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Set t ■ At, n • 1. 


n 


(2) The fluxes are co^nited using u.. . 

(3) Tlje fluxes are updated on those parts of the boundary where is 
specified. 

(4) The fluxes are differentiated to obtain the "parabolic" teim 
3x 3y* 

(5) is computed from equation (6). 

(6) is updated on those parts of the boundary where u is 
specified. 

(7) t t ♦ At , n ■*- n ♦ 1 , go to step (2) . 


Note that only first-order (numerical) derivatives are required in 
steps (2) and (4). Thus there is one derivation routine, which is 
the usual straightforward T^ formula. In certain cases, however, 
high mode Interactions cause instabilities, and the high coefficients 
generated by differentiation must be smootKed. The smoothing: 
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is then included in the derivative routine. The updating of values in 
step (3) varies from problem to problem. The best results have been 
obtained by overspecifying data, i.e. prescribing f when is' given. 

It is also possible to compute f, g from the correct data while 

the tangential derivative needed is extrapolated along characteristics. 

Of course, the characteristics pertain to the wave equation [eq. (5)]. 


The analysis of reference 1 could be applied directly to equation (3) 
to prove stability, except for the very complicated structure of the 
"parabolic term" operator, which also depends on boundary conditions. 

For certain cases, Orszag (ref. 2) has numerically computed the 
eigenvalues of this operator, which turned out to be all real — a 
sufficient condition for stability, according to Gottlieb and Gustafsson 
(ref. 1). 


4 



Actual conputation, as presented in the following sections, amply 
justifies the assumption that generalized Dufort-Frankel Tschebyscheff 
schemes are stable for appropriate multipliers y. 

DIFFUSION EQUATIONS WITH CONSTANT COEFFICIENTS 


This section suanuirizes results obtained for the equation 


u u ♦ 2eu 
XX yy xy 


( 11 ) 


For parabolicity, |e| < 1. The dcmiain considered is |x 111 . |y| 11 . 
t 0, with u specified initially and either u or specified 
at |x| » 1, |y| » 1. 


An analytic solution has always been used to check accuracy, and also 
to supply initial and boundary values, as well as first level values u(x,y,t » 
at) . No smoothing of coefficients was done, and fluxes were overspecified 
on boundaries. The first check was to test whether our numerical method 
[eq. (10)] actually admits larger At than the obvious explicit 
method : 
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A comparison, done for e = 0 and Oivichlet conditions everywhere, 
showed that at may be easily taken to be up to 50 times larger than 
the stability limit of equation (12). Thus, equation (10) with y ~ S 
converges for at * 0.001, N * 32, while the explicit method diverges 
at N « 32, at * 0.00005 and converges only at N * 32, At = 0.00002. 

Next the various values for the multiplier y were considered. 

In most of these computations, y 4 is the stability limit. For 
instance, N > 32, At « 0.001 will diverge at y 3, but will converge 
at Y * 5; N = 16, At * 0.001 will converge at y = 3.375 but diverge 
at Y • 3.25. Of course, the smaller y. the better the approximation; 
however, the error increases only slightly with y> 

Another test is the sensitivity of the method to the definition of 

ax., ly.. If all these quantities are taken as equal [as opposed to eq. (4)], 

) ^ 

the results are less accurate and the multiplier y must be larger to 
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produce stability. On the other hand, the ambiguity in flux definition 
has no influence — even if fluxes must be updated on the boundary. 

The boundary flux treatment is evidently the most complex of the 
questions connected with the present scheme. Overspecification produces 
satisfactory results. The characteristic extrapolation, i.e., 

u^ + eUy » u^ at X = ±1 (13) 

seems to be less accurate. Here an analytic result would be extremely 
useful . 


CONNECTION-DIFFUSION PROBLEMS 

In this section the following two cases are considered: 

(1) Constant coefficient problem: 

a ua va > ^ (a * a l 

^y Pe '■’^xx 

0 £ x,y £ IT, t ^ 0 

• a given on y * 0,ir 
a = 0 on X * 0,v 

a given at t = 0 (14) 

(-) Temperature distribution in a Benard cell: 

♦t • U(x,y),^ . V(x,y),y . i 

0 < x,y < u, t ^ 0 
a*0ony*0, a = l ony=* 

a^ = 0 on X = 0,1^ 

a(x,y,t = 0) given (15) 

The speeds U(x,y), V(x,y) representing a vortical flow are prescribed: 

U(x,y) = -cos y sin x 

V(x,y) » sin y cos x (16) 
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The quantity Pe is the Peclet nximber of the problem: Pe • 


typical speed x typical length 
conductivity 

specified as 


f 


Pe ’x 




In both cases j the fluxes have been 


g 





which is appropriate even for variable coefficients, since (U.V) is 
divergence free. The flux fg^ = - !)♦ was prescribed on the boundary 
carrying Neumann conditions. Initial conditions have been generated 
arbitrarily for both equations (14) and (15), and the convergence to a 
steady solution was tested. 


( 17 ) 


For the constant-coefficients case, an analytic solution is easily 
obtained which is used as before to supply boundary conditions and verify 
accuracy. For instance, one may take 


» (cos X - 


UPe 


sin x) sinh 


W 


(U^ » y2)Pe 


- *)• 


Kf * f ) 


(18) 


as a steady solution of equation (14). Formula (18) clearly shows the rapid 
variation of 4 with y as Pe increases. Indeed, at U*V»1, Pe » 5, 
the computational results became quite poor, since the few modes used 
(To, T; . . . Tg or eventually Tq, Ti . . .Tig) could not properly 
resolve = e^*®X. For smaller Pe (larger conductivity), very accurate 
results were readily produced. 


For equation (15) collocation was used to deal with the variable 
coefficients. It appears that enough high modes were generated to distort 
the solution. Some plots even show minima and maxima (there should be 
none inside the domain) at the nodal points. The smoothing formula [eq. (10)] 
provides a simple remedy to this unwanted phenomenon. For Pe = 10, N = 16 
successful computations have been performed, with results in complete 
agreement with numerical data obtained otherwise. (There are no analytic 
solutions for this case). Note that the Peclet number seems higher, but, 
in fact, because U, V vary inside the domain, the boundary layers of 
equation (15) are less abrupt than those of equation (14) and resolution 
is easier. Using polynomials of higher degree (say 64), one could treat 
even higher Peclet numbers. (See also figures 1 and 2). 
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The convergence to the steady state is quite rapid, as seen in figures 
3 and 4. Analytic solutions of equation (14), say, with homogenous conditions, 
decay at least as fast as 




(U^ y2)pe 


Thus, one may use time-dependent equations — or even artificially time dependent 
equations — to solve steady-state convection-diffusion problems efficiently. 


BURGERS EQUATION 


This is an experiment in treating nonlinear equations. 




< 


to be solved is 


u + uu 
t s 


3^u 
as 2 


where s is defined by 
ooc ♦ 8y 

s a • i.'.l s-. o, 8 constants. 


(19) 


( 20 ) 


Of course, equation (19) is rewritten as an equation in (x,y,t) domain, 
with boundary and initial conditions. However, these are specified in 
such a way as to agree with equation (19). Specifically, we solve 


“t * ' '’^“^“xx * ®^”yy * 


-1 1 x,y £ 1 , t ^ 0 

u(x,y,t « 0)» Ug(s) 

^[boundary * [boundary’ 

The fluxes are defined as 


f » v(a^u^ * aSUy) - y 

g » v(a8Uj^ ♦ 8^Uy) - I u2 (22) 

This problem has been selected because of its explicit solutions, 
which enable comparison with the numerical results. Satisfactory convergence 
is again easily obtained. An example is presented in figure 5. An 
interesting point is that equation (21), although nonlinear, shows better 
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results when no smoothing is applied to the high coefficients (cf . ref. 4) . 
This is probably due to the fact that there is enough viscosity in the 
equation itself. 



which is a solution of equation (19) with a » 0.7, 6 « 0.8, v « 0.1. 

At time t * 0, the maximum value of u is ^ 3 while the average value of u 
over -1 £ x,y ^ 1 is « 0.2. Thus there are steep boundary layers to be 
resolved. If the smoothing (eq. (10)] is applied, the result is much too 
danq>ed; without smoothing, for N « 16, y • 15, an error of about S • 10'*^ 
is maintained for 0 < t < 0.5, corresponding to 500 steps (At « 0.001). 


DISCUSSION 


The purpose of this research program was to obtain better time 
increments than those allowed by the naive scheme [eq. (12)], while 
retaining the desirable spectral accuracy. Here other possible solutions 
to this problem will be briefly described. One is presented in reference 5, 
where the explicit time-stepping operator is modified in such a way as to 
make it bounded for all At. This method is efficient for Fourier expansions, 
but not for Tschebyscheff expansions. After some experimentation it was 
found that the modified operator, while stable, differs too much from the 
exact one, and the approximations are inconsistent. 

Another method is to employ approximate inverses (ref. 2) . These 
are operators, simple in structure, which approach the spectral ones. 

Implicit time stepping is used, with arbitrary At, and the inversions 
needed are performed on the simpler operators. To make the procedure 
efficient in multidimensions, some type of splitting must be used, introducing 
problematic boundary treatment. Mixed terms also pose problems (ref. 6) 
coBQ)ared with the simplicity of equation (8). Yet another possibility is 
to generalize Saulev’s scheme (ref. 7) to spectral methods in the same 
way equation (2) is generalized to equation (3). This seems a very promising 
perspective, and a probable direction for future work. 
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CONCLUSIONS 


A siaple and efficient time-advancing scheme has been presented 
for spectrally solving parabolic equations with general boundary conditions. 
Various equations have been used, obtaining satisfactory results, especially 
for convergence to steady states and Divichlet boundary conditions. 

While the main theoretical question ~ stability — is still open, 
these preliminary results provide a fine foundation for further research. 
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Figure 1. Divergence of B^nard problem at Pe = 50 showing lines of constant ♦ at 
t = 4, computed using nine modes (Tq to Te), y = 5, At = 0.01. 
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Figure 2. Divergence of Benard problem: exact steady solution obtained from a 

multigrid program. (The x and y axes are interchanged with respect 
to figure 1.) 
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Figure 3 


Convergence to steady state: solution of equation (14) showing arbitrary 

initial data; nine modes used (Tq to Ta) , y= 3. At =0,01, U=V= 1, 

Pe * 1. 
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Figure 4 


Convergence to steady state: solution of equation (14) with arbitrary 
initial data showing the function at t > 2; nine modes used (Tq to Te) , 
Y = 3, At = 0.01, U = V = 1, Pe * 1, maximum error < 1 x 10~**. 
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Figure 5. Solution to Burgers equation at t = 0,5 with error 5 x 10“^; 1/ modes 
used (Tq to Tjg), Y = 5, At = 0.001, a = 1, 6 = 2, v = 0.001. 

Note that although equation (21) is written in terms of x and y, the 
solution u = constant on lines s = constant (straight lines of slope 
- 1 : 2 ). 
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